;; LTDL_LIBRARY_PATH=/opt/local/lib
;; GUILE_LOAD_PATH=.

;;(import (srfi srfi-1))

(import (rnrs)
        (grid-gen common)
)

;; (module-on real-8)
;; (module-on OpenCL)
;; (module-on OpenCL-CPU)

(module-on fortran77)
(module-on HD)
        
;; (module-on Lax-Friedrichs)
(module-on Roe)

(module-on Einfeldt)

;; (module-on Harten2)

(module-on Osher)

(module-on grid)
(module-on maxima)
(module-on OpenDX)
        
(module-on OpenMP)
        
;; (module-on MPI 8)

;; (Harten2:set-eps! 1.0e-6)

(def-float tim 0.0)
(def-float tau)

(define grid-dims '(200 400))
;;(define grid-dims '(1000 2000))

(define pi 3.1415926)

(def-float xmin -1.0)
(def-float xmax  1.0)

(def-float ymin  0.0)
(def-float ymax  4.0)

(def-float dphi  0.25)

(def-float phimin -0.25)
(def-float phimax  0.25)

(def-float rmin 1.0)
(def-float rmax 5.0)

(def-float gkf  4.0)

(def-float tmax 0.8)

(def-float sr 0.1)
(def-float sl 1.0)

(def-float rho0  1.0)
(def-float p0    1.0)
(def-float u0    0.0)
(def-float v0    0.0)

(def-float rho1   10.0)
(def-float p1     1.0)
(def-float u1     0.0)
(def-float v1     5.0)

(def-float x0     0.0)
(def-float y0     2.0)
(def-float r0     0.1)

(def-float cour   0.1)

(define gr (grid:cart grid-dims))

(define gr-cond `(< (abs ,(gr 'x)) ,sr))
;;; (define gr-cond `(< (abs (- ,(gr 'y) 2.0)) ,sr))

(define bul-cond
       `(< (sqrt (+ (** (- ,(gr 'x) ,x0) 2.0)
                    (** (- ,(gr 'y) ,y0) 2.0)
                 )
           )
          ,r0
        )
)

(gr 'add-bnd 'y-
   (beg
     (set-val (gr 'r)   `(if ,gr-cond ,rho1 ,rho0))
     (set-val (gr 'p)   `(if ,gr-cond ,p1   ,p0))
     (set-val (gr 'u)   `(if ,gr-cond ,u1   ,u0))
     (set-val (gr 'v)   `(if ,gr-cond ,v1   ,v0))
    )
)

(gr 'add-corr 'smooth)
(gr 'add-corr 'err)

(define grdf1 (lambda (ct cr z)
  (let ((phi `(* 2.0 ,pi ,dphi (- ,ct 0.5)))
        (r   `(+ ,rmin (* (- ,rmax ,rmin) ,cr)))
       )
    (list `(* ,r (sin ,phi))
          `(* ,r (cos ,phi))
          z
    )
  )
))

(define grdf0 (lambda (x y z)
  (list `(+ ,xmin (* (- ,xmax ,xmin) ,x))
        `(+ ,ymin (* (- ,ymax ,ymin) ,y))
        z
  )
))

(define grdf01 (lambda (x y z)
  (list `(+ ,xmin (* (- ,xmax ,xmin) 0.25 (- (** 5.0 ,x) 1.0)))
        `(+ ,ymin (* (- ,ymax ,ymin) 0.25 (- (** 5.0 ,y) 1.0)))
        z
  )
))

(define grdf3 (lambda (x y z)
  (list `(+ ,xmin 
             (* (- ,xmax ,xmin) 
               ,x
                (+ 1.0 (* 4.0 (- 1.0 (/ (- ,y ,ymin) (- ,ymax ,ymin)))))
             )
         )
        `(+ ,ymin (* (- ,ymax ,ymin) ,y))
         z
  )
))

(define grdf2 (lambda (x y z)
  (list `(+ ,xmin
            (* (- ,xmax ,xmin) ,x)
            (* ,gkf (* (* 2.0 (** ,y 2.0)) (- 1.5 ,y)) (- ,x 0.5))
         )
        `(+ ,ymin (* (- ,ymax ,ymin) ,y))
         z
  )
))

(define shc (lambda (v1 v2)
  `(if (and (< (abs ,(gr 'x)) ,sr) (< (- ,(gr 'y) ,rmin) ,sl)) ,v1 ,v2)
))

(define shc1 (lambda (v1 v2)
  v2
))

(def-float G 3.0)

(define mypot (lambda (x y z)
  `(* ,x ,G)
))

;;(set! PHYS:Potential mypot)

(gr 'set-geom grdf0)

(code-gen "stream"
  (beg
    (gr 'init-grid)
    (gr 'init (shc1 rho1 rho0) (shc1 p1 p0) (shc1 u1 u0) (shc1 v1 v0))

;    (gr 'init (if bul-cond rho1 rho0)
;              (if bul-cond p1 p0)
;              (if bul-cond u1 u0)
;              (if bul-cond v1 v0)
;    )
    (do-while (and (< tim tmax) (= (gr 'err) 0))

      (gr 'update-bnd)
      
      (if (= (gr 'err) 0)
        (beg
          (gr 'CFL cour)
          (if (> (+ (gr 'tau) tim) tmax)
              (set-val (gr 'tau) (- tmax tim))
          )
        )
      )
      (main-p
        (do-every 1 'sec
          (print "t=" tim "tau=" (gr 'tau))
        )
      )
      (gr 'calc)
      (if (= (gr 'err) 0)
          (set-val tim (+ tim (gr 'tau)))
          (beg
            (set-val (gr 'tau) (* 0.5 (gr 'tau)))
            (print "new Tau=" (gr 'tau))
          )
      )
    )
    (OpenDX:save gr "stream.dx")
  )
)

(display "\n")
